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ABSTRACT 

We study the long term evolution of magnetic fields generated by a collisionlcss relativistic e + e~ shock 
which is initially unmagnetized. Our 2D particle-in-cell numerical simulations show that downstream 
of such a Weibel-mediated shock, particle distributions are close to isotropic, relativistic Maxwellians, 
and the magnetic turbulence is highly intermittent spatially. The non-propagating magnetic fields 
in the turbulence form relatively isolated regions with transverse dimension ~ 10 — 20 skin depths. 
These structures decay in amplitude, with little sign of downstream merging. The fields start with 
magnetic energy density ~ (0.1 — 0.2) of the upstream kinetic energy within the shock transition, but 
rapid downstream decay drives the fields to much smaller values, below 10~ 3 of equipartition after 
~ 10 3 skin depths. 

In an attempt to construct a theory that follows field decay to these smaller values, we explore the 
hypothesis that the observed damping is a variant of Landau damping in an unmagnetized plasma. 
The model is based on the small value of the downstream magnetic energy density, which suggests that 
particle orbits are only weakly perturbed from straight line motion, if the turbulence is homogeneous. 
Using linear kinetic theory applied to electromagnetic fields in an isotropic, relativistic Maxwcllian 
plasma, we find a simple analytic form for the damping rates, 7^, in two and three dimensions for 
small amplitude, subluminous electromagnetic fields. We find that magnetic energy does damp due 
to phase mixing of current carrying particles as {uj p t)~ q with q ~ 1. This overall decay compares well 
to that found in simulations, since it depends primarily on the longest wavelength modes, kc/u> p <C 1. 
However, the theory predicts overly rapid damping of short wavelength modes. We speculate that 
magnetic trapping of a substantial fraction of the particles within the highly spatially intermittent 
downstream magnetic structures may be the origin of this discrepancy. In addition, trapping may 
form the basis for MHD-like behavior, permitting a small fraction of the initial magnetic energy to 
persist for times much greater than have been followed in the simulations. 

We briefly speculate on other physical processes, which depend on the presence of suprathermal 
particles, that may cause the generation of longer wavelength magnetic fields that create a magnetized 
plasma (krLarmor <C 1), in which the damping is not as fast. However, absent such additional physical 
processes, we conclude that initially unmagnetized relativistic shocks in electron-positron plasmas are 
unable to form persistent downstream magnetic fields. These results put interesting constraints on 
synchrotron models for the prompt and afterglow emission from GRBs. We also comment on the 
relevance of these results for relativistic shocks in electron-ion plasmas. 
Subject headings: shock waves - turbulence - gamma ray: bursts - plasmas 
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1. INTRODUCTION 

The prompt emission and afterglows of gamma-ray 
bursts (GRBs) may be manifestations of ultrarelativistic 
shock waves, propagating in media where the large scale 
upstream magnetization is too weak to affect the shock 
structure, and too weak, if simply compressed by the 
shock, to provide the magnetization inferred from syn- 
chrotron models of the burst emission (see Piran 2005ab 
and references therein). The weakly magnetized outflows 
in the rotational equators of rotation powered pulsars 
(Coroniti 1990) may also be sites of essentially unmagne- 
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tized shock waves terminating the relativistic winds in a 
region which occupies a finite latitude band with respect 
to the rotational equator of the underlying neutron star. 
Yet the radiation from these systems has been modeled 
as being due to synchrotron radiation, which requires 
the presence of a magnetic field strong enough to deflect 
particles through a substantial fraction of their Larmor 
orbits. At the very least, this requires magnetic struc- 
tures with amplitude, SB, whose characteristic dimen- 
sions, Rb, are comparable to or larger than mc 2 /e8B, 
for all particle energies inferred to contribute to the ob- 
served radiation. 

Unmagnetized anisotropic plasmas spontaneously gen- 
erate small-scale magnetic fields via the Weibel instabil- 
ity (Weibel 1959). Shocks have strong plasma anisotropy 
in the transition layer separating the upstream and 
downstream media, as well as in the foreshock where 
downstream particles escape into the upstream, which 
provides the free energy for generating magnetic fields 
with spatial scales on the order of the plasma skin depth, 
c/oj p , where uj p is the plasma frequency. Medvedev 
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and Loeb (1999) and Gruzinov and Waxman (1999) ar- 
gued that the relativistic form of this instability (Yoon 
&; Davidson 1987) could macroscopically magnetize the 
downsteam plasma where the GRB radiation arises (Lar- 
mor radii small compared to flow scale). Through numer- 
ical and analytic studies, they and various authors (Silva 
et al. 2003, Medvedev et al. 2005) have argued that the 
instability forms filaments of electric current and B field. 
These filaments merge and cause magnetic energy to cas- 
cade from the initial microscopic scale ~ c/lu p to larger 
scales. Thus, the filamented plasma becomes magnetized 
with a B field hypothesized to survive throughout a large 
fraction of the shocked medium. 

Such inverse cascades have been observed in 2D and 
3D particle-in-cell simulations in the foreshock region, 
the part of the shock structure where the upstream and 
downstream media interpenetrate and the stream fila- 
mentation modes grow from noise with most incoming 
particles still undeflected from the upstream flow. Sim- 
ulations show that current filaments merge and grow in 
amplitude until they reach the magnetic trapping limit, 
where the filament currents and their magnetic fields 
become comparable to the Alfven limit (Davidson et 
al. 1972; Kato 2005; also see Milosavljevic, Nakar, & 
Spitkovsky 2006; Milosavljevic & Nakar 2006a). At that 
point the particles' orbits on the filament boundaries be- 
come chaotic, the filaments disorganize, and scattering 
from the disorganized magnetic fluctuations halts the 
streaming of the bulk of the plasma, isotropizing and 
thermalizing the flow, all within a layer tens to hundreds 
of skin depths thick (Spitkovsky 2005), in accord with 
Kato (2005) 's model for shock formation. The magnetic 
energy is as high as 10-20% of the bulk plasma flow en- 
ergy within this scattering layer, where the density jump 
between upstream and downstream occurs. 

These initially small-scale B-fields must survive for 
tens of thousands to millions of inverse plasma peri- 
ods to serve as the source of the magnetization invoked 
in synchrotron models of GRB emission (Gruzinov & 
Waxman 1999; Piran 2005ab; Katz, Keshet, & Waxman 
2007). Long-lived fields are also required in Diffusive 
Fermi Acceleration (DFA) models used to explain the 
appearance of the nonthermal particle spectra observed 
through their synchrotron and inverse Compton emis- 
sion in GRBs, and in pulsar wind nebulae (PWNs) and 
other sites of nonthermal photon emission in relativis- 
tic flows. However, present simulations (Kazimura et al. 
1998; Silva et al. 2003; Frederiksen et al. 2004; Hededal 
et al. 2005; Nishikawa et al. 2003, 2005) have not fol- 
lowed the flow through the shock transition layer into the 
downstream region because they employ periodic bound- 
ary conditions, are too small, or both. The question of 
the structure and long-term survival of the B-fields has 
remained open (see for instance, Gruzinov & Waxman 
1999; Gruzinov 2001ab; Medvedev et al. 2005). 

In this paper, we show via linear theory that the 
phase mixing between individual particles and organized 
macroscopic currents implies rapid decay of the mag- 
netic energy in the downstream medium. We begin by 
briefly describing the essential features of the plasma 
produced by the shock via numerical simulations in SjSJ 
These simulations show the downstream magnetic struc- 
tures are non-propagating in the frame of the down- 
stream medium, and are intermittent spatially, organized 



into clumps and filaments of magnetic field with typi- 
cal diameter ~ 10 — 20 skin depths, immersed into a 
highly isotropic plasma. In this downstream isotropic 
Maxwellian particle distribution, we calculate the damp- 
ing rates from Vlasov linear response theory, assum- 
ing the particles are unmagnetized (Larmor radii S> 
magnetic clump diameter), recovering a result due to 
Mikhailovskii (1979) in $3]and Appendix A. Taking snap- 
shots from the simulations as initial conditions, we cal- 
culate the decay of the magnetic field as a function of 
time and position, and compare these theoretical calcu- 
lations with simulations. We find that the theory does 
reasonably well in estimating the decay rate of the total 
magnetic energy as t~ q with q ~ 1. However, it overes- 
timates the damping rate of shorter wavelength modes. 
We speculate in ^4]that this discrepancy may result from 
the magnetic trapping of a large fraction of the particles, 
which suppresses the disorganizing effect of phase mixing 
on the currents. We also discuss the effects of an inverse 
cascade on the persistence of magnetic fields, though we 
find little evidence for its relevance from numerical simu- 
lations. We summarize our results and draw some impli- 
cations for Fermi acceleration and for the magnetization 
required for synchrotron emission in models of GRBs and 
other systems in $5] 

2. SIMULATION RESULTS 

Spitkovsky (2005) and Spitkovsky and Arons (in prep) 
describe a series of 2D and 3D simulations of relativistic 
shock waves in e + e~ plasmas, for varying values of the 
upstream magnetic field B\, including B\ — 0. These are 
Particle- in-Cell (PIC) simulations (Birdsall and Langdon 
1991), using the code TRISTAN-MP, developed by one 
of us (AS). It is a heavily modified descendant of the 
publicly available code TRISTAN (Buneman 1993). We 
simulate shocks by injecting cold relativistic plasma par- 
ticles at one end of a large domain and allowing the par- 
ticles to reflect off a fixed conducting wall at the other 
end of the box. In this paper, we present 2D calculations 
utilizing boxes as large as 50,000 x 2048 cells with up 
to 1.35 x 10 10 particles which allows us to fully resolve 
shock formation. Although we track all three velocity 
and field components, due to the two-dimensional sym- 
metry only particle velocities in the plane of the simu- 
lation are non-zero for initially cold flow, and only the 
out-of-plane component of the magnetic field is excited 
by in-plane currents. In-plane electrostatic fields are also 
included. 

The interaction of the reflected pair plasma with the 
incoming stream forms a shock, which propagates toward 
the plasma injection surface. In the simulations, one 
plasma skin depth spans 10 cells based on the upstream 
parameters, i.e., X s i = c/u) p x = y/ m±c 2 7i/47re 2 ni, 
where n\ is the upstream total density of electrons and 
positrons (as measured in the frame of the simulation), 
ojpi is the upstream plasma frequency, and r )\m±c 2 is 
the upstream flow energy/particle. The time step is 
At = l/20tL> p i. In the upstream skin depth units, the 
largest boxes are then 5000 x 205 c/u> p i. The longest 
simulation was evolved for 5300^"^. We have checked, 
by using boxes of increasing transverse dimension, that 
the periodic boundary conditions used on the transverse 
coordinate do not affect the scale of the magnetic struc- 
tures formed. In the results exhibited below, we use 64 
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particles per cell (32 per species) in order to suppress 
particle noise, which enables us to follow the dynamics 
of small amplitude fields. We have performed conver- 
gence studies, varying the number N of particles per cell 
from 4 to 64, confirming that the noise level decreases as 
1 / V~N, while the gross qualities of the shock structure re- 
main the same. Figure [T] shows the snapshots of density 
and magnetic energy from a typical 2D simulation. 7 Co- 
ordinates are in units of the upstream skin depth, cjw p \. 
In the simulation shown, the upstream flow moves to the 
left with 7i = 15. 

Our simulations are large enough to permit the com- 
plete development of the shock and show the main 
features of contemporary collisionless shock simulations 
even in reduced dimensionality. For instance, we see 
the factor of ~ 3.13 increase in density between the up- 
stream and the downstream (Fig. [TJ:), which is the ex- 
pected compression factor for this 2D plasma when the 
plasma properties are measured in the rest frame of the 
downstream plasma (Gallant et al. 1992; Spitkovsky and 
Arons, in prep). Current filaments (geometrically, these 
are out-of-plane sheets) show up as the enhancement in 
plasma density and magnetic energy density in the fore- 
shock (Fig. []Jb). The scale of the filaments grows to- 
wards the shock through merging. The shock is located 
where the density filaments completely merge and are re- 
placed by a quasi-homogeneous medium. The subject of 
the saturation of the Weibel instability will be explored 
in greater detail by Spitkovsky and Arons (in prep). 

At the shock transition layer, the filaments disorga- 
nize and become clumps of magnetic energy (in 2D the 
only appreciable magnetic component is the out of plane 
B z ). Note that these magnetic clumps lose intensity 
the further downstream they are from the shock (Fig- 
ure [T]d) . The magnetic energy peaks before the density 
completes its rise (as we see in comparing c and d in 
Fig. Q}, i.e., the instability saturates at the Alfven criti- 
cal current before the shock fully develops. We also note 
that the particle distribution function changes from an 
anisotropic free-streaming population to an isotropic (in 
the downstream rest frame) thermal population. In the 
downstream medium, we find that the difference between 
the perpendicular and parallel momentum is extremely 
small, < 1%. 

As Figures Q]b an d d suggest, even though locally in 
the shock magnetic energy can reach close to equiparti- 
tion (epsiloriB ~ 15% when averaged over the transverse 
dimension) , the magnetic fields decay in the downstream 
region of the shock. To study this in further detail, we 
present several cross sections of the plasma at times sep- 
arated by 450 ujpi in Figure O Since our simulation is 
performed in the downstream frame (frame of the reflect- 
ing wall) , we see the shock propagate through the box at 
rs 0.5c (value appropriate for 2D relativistic gas). 

Downstream from the shock, the magnetic clumps 
weaken as a function of time and appear to be non- 
propagating. In Figure [2 we see this in panels a-d and 
in the black through blue curves in panel e, which show 
magnetic energy averaged over the transverse dimension 

7 Note that we have defined the magnetic energy fraction eg = 
B 2 /47T7i ni mc 2 in terms of the upstream kinetic energy density, 
instead of the downstream thermal energy density as is a common 
practice in the GRB community (c.f., Piran 1999). 



of the simulation at times from panels a-d. Although the 
separation between prominent clumps that have larger 
field strength seems to increase with time, this is due to 
the faster disappearance of small-scale clumps, presum- 
ably due to decay, while the location or size of strong 
clumps does not significantly evolve. There is little evi- 
dence for clump merging far from the shock. In Figure 
[2] we also plot a vertical line showing the location where 
we will decompose the magnetic fields into Fourier modes 
later in §3, in order to study the comparison of the theory 
of B-field decay with the numerical experiments. 

Finally, we mention that 3D shock simulations also 
show similar behavior (including the lack of substan- 
tial downstream merging of magnetic structures) for the 
shorter times that can be followed with present 3D sim- 
ulations. Topologically, the magnetic clumps in 2D be- 
come large looping structures in 3D. If our 2D plane is 
thought of as a slice through a 3D simulation, the 3D 
loops connect field emerging from one clump and return- 
ing to another. The orientation of 3D loops would be 
mostly perpendicular to the original direction of the flow. 
The particle distribution function is also an isotropic, rel- 
ativistic Maxwellian in the downstream region. 

3. DOWNSTREAM EVOLUTION OF MAGNETIC 
TURBULENCE 

The simulations summarized in <J2] suggest that mag- 
netic turbulence decays in the isotropic post-shock 
plasma. We now attempt to understand this theoreti- 
cally with the goal of finding a model that allows ex- 
trapolation of the magnetic evolution beyond the length 
of time that can be studied in direct shock simulations. 
The simulations show that the downstream plasma is 
isotropic and the downstream particle distribution func- 
tion is well described by a relativistic Maxwellian. For 
the purposes of this section, we assume the downstream 
field amplitudes are so small that magnetic trapping is 
unimportant for almost all particles; therefore, their or- 
bits are almost straight lines. We will revisit this as- 
sumption in 2) 

We begin by deriving the linear plasma response for 
electromagnetic fluctuations in an isotropic relativistic 
plasma with small field fluctuations (Mikhailovskii 1979). 
The Vlasov equation for each species is 



at m e V c 



V p / S = 0. (1) 



Here, s is the species label, + for positrons and — for 
electrons, with q± = ±e. We linearize this equation with 
fs — > fos + Sfs, B — > SB and E — > SE, with the initial 
conditions downstream of the shock that tell us that SE 
and 5B are small in the sense that (SE 2 + SB 2 ) /8nnT <C 
1. Then Sf s /fos is also small. The linearized Vlasov 
equation is 

^ +v VSf s + ^(se + -xSb)- V p f Qs = 0. (2) 

at m e V c / 

The plasma couples to the field through the Maxwell 
equations 

1 dSB 



V X SE: 



V xSB -. 



c dt ' 

4ir sj + 1 dSE 

c c dt 



(3) 
(4) 
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Fig. 1. — Snapshot of a region from a large 2D relativistic shock simulation. Incoming 7 = 15 flow is moving to the left, while the shock 
moves to the right at ?s 0.5c. The postshock plasma is stationary in the frame of the simulation, a) Density structure in the simulation 
plane showing the plasma density enhancements in the foreshock region that steadily grow up to the shock transition region, where the 
density becomes homogeneous. Density is normalized in units of the plasma density far upstream. The density jump, (n2)/{ni) = 3.13, 
shown is exactly what is predicted by the hydrodynamic jump conditions for 71 = 15 in a 2D gas. b) Magnetic energy, normalized in 
terms of upstream energy of the incoming flow: eg = B 2 /A-K^inimc 2 . The upstream magnetic filaments, which can be visualized as sheets 
coming out of the page, that are formed by the Weibel instability reach a peak just before the shock. These filaments become clumps of 
magnetic field, which can be visualized as cross-sections of loops that are transverse to the page in the downstream region, where they 
slowly decay away. Note that the B-field, which is organized into upstream filaments or downstream clumps, always points in or out of 

the page. A power law scaling, e^/ 4 , was applied to stretch the color table to show weak field regions; this is reflected in the colorbar. c) 
Plasma density averaged in the transverse direction as a function of the distance along the flow, d) Magnetic energy density averaged in 
the transverse direction, as a function of distance along the flow. 



where Sj = J2 S ? s / v ^fsd 3 p is the current density. 

We orient coordinates so that 5E is along x, the wave 
vector lies along y, and SB is along z. Assuming a wave 
solution with amplitudes oc exp (— iu>t + iky), equation 
© becomes 



i (kv v - uj) Sf s + — (8E + ^-Sb) 
m P \ c J 



dfos Qs v x dfos 



dp x m e c dp y 
(5) 

Using equation ^ to eliminate SE yields 



Sfs 



. q s 6B ( dfo s kv x df 0s 



(6) 



m e kc \ dp x lu — kv y dp y 
Using (j4]) yields 

i (uj 2 - k 2 c 2 ) 5B = -iirkcSj. (7) 
Using ([6]) in the definition of 5j yields 

Sj = -i X ujSE (8) 



where \ is the plasma susceptibility, 



4tt X = 



p,NR,s 
W 2 



_d_ 

dp x 



kv T 



LU — kVy dp y 



(TP. 

(9) 

Here n s is the number density of the electrons (s = — ) or 
0- positrons (s = +) in the downstream plasma, w Pi nr.s = 
\J 4irq2n s /m e is the non-relativistic plasma frequency. 
We will assume by charge neutrality that in the down- 
stream plasma /o+ = /o- and n + = rt_ so that the sum 
over equal mass species is trivial. 

The dispersion relation for normal modes in the plasma 
plus electromagnetic field is 



2^,2 



k 2 c 



(1 + 4ttx)=0. 



(10) 



However, we emphasize that the linear relation between 
the currents and the electric field (eq.[8]) through the 
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Fig. 2. — Cross-sectional views of the magnetic energy density for the same region as Figure[T] but at different times separated by 450 
ujp 1 (a-d). The time of Figure [T] corresponds to panel b) in this plot. In panel e) we plot the B-field energy, averaged over y-axis as a 
function of distance along the flow direction with one curve for each of the top panels: (a) black, (b) red, (c) green, (d) blue. We also plot 
a vertical line at x = 840c/ u> p to define the region over which we perform a modal analysis in i[3land Figure [3] 



susceptibility (J9j) applies to all (small amplitude) fluctu- 
ations, not only to normal modes. 

We evaluate equation © for distribution functions 
that are isotropic in two and three dimensions in Ap- 
pendix A. For a two-dimensional isotropic distribution 
fu nction, we write f 0s (p) = J{p2d)g{Pz) 1 where p 2 d = 

\jv\ + Py an d perform the integral over p z such that 

J g(p z )dp z = 1. For a three-dimensional isotropic dis- 
tribution function, we write fo s (p) = fos(p3d), where 



P3d = yPx+Py+Pz- I n the three-dimensional case, 

we confirm Mikhailovskii's (1979) result for a relativis- 
tic isotropic plasma. As we show explicitly in equation 
(|A5[) . subluminal waves (u> r < kc) are damped, where 
uj r = 0?(w). The basic physics is that of Landau damp- 
ing (Stix 1992) or more precisely, phase mixing. For a 
subluminous wave in an isotropic relativistic plasma with 
oj r /k < c, there always exists some angle <f) between the 
wave vector k and a particle's momentum where a par- 
ticle moving with speed c(3 is in phase with the field 
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component, i.e., /3cos0 = u> r /kc. 

When the phase velocity of the field fluctuations is very 
small compared to the thermal speeds of the particles 
(which includes the zero wave velocity case, SR(o;) = 0), 
Landau damping takes on the character of simple phase- 
mixing. The currents, which support the field fluctu- 
ations, are composed of particles. These current car- 
rying particles have random motions which carry them 
out of the magnetic structure, which the currents ini- 
tially support. As a result, these currents and their mag- 
netic field fluctuations are disorganized (or damped) on 
the transit times of the particles across the field fluc- 
tations. It should be emphasized that the formal the- 
ory includes both this phase mixing limit of the damp- 
ing process and the limit in which the phase 4-speeds 

(3 P / ^Jl — (3p, (/3 P = 9?cj/fcc), of the fields' Fourier com- 
ponents are high compared to the particles' random 4- 
speeds (3/ y/ 1 — [3 2 . The applicable limit depends on 
whether the wave phase 4- velocity is large or small com- 
pared to the mean 4-velocity of the particles in the dis- 
tribution function. Hammett, Dorland & Perkins (1992) 
present a simple picture of this process in the case of 
non-relativistic electrostatic waves. Arons, Norman & 
Max (1977) show how phase mixing works for electro- 
magnetic waves in a relativistic Maxwellian plasma. In 
any case, the end result is that there is a net power flow 
from fields to particles, i.e. the fields decay, which we 
will demonstrate explicitly below. 

We now study the action of such a plasma on an ar- 
bitary field of initial B-field pertubations as generated 
by the Weibel mediated shock. We evaluate equation ((9]) 
numerically for two and three dimensions setting oj r = 0, 
because of the non-propagating nature of the magnetic 
clumps, which we infer from visual inspection of the sim- 
ulations. In the limit k <C ^p/c, we find: 



4ttx 



'4 \kj, 



2D 
3D 



(11) 



Note that the 2D and 3D results only vary by a numerical 
factor. Hence, long wavelength modes have the same 
qualitative behavior in two and three dimensions. 

The plasma susceptibility, which we calculate in Ap- 
pendix A (see also eq.pj]), can be utilized to calculate 
the evolution, i.e., the damping or growth, of an initial 
field of fluctuations. Appendix B contains this calcula- 
tion in more detail; the result is 



d\6B k 
dt 



-2-Yk\6B k \ 



(12) 



where 7^ = (kc) uj~ 3(47rx)~ (see also eq. |B8| ). The 
asymptotic forms of x from equation (fTTj) for 3D and 2D 
gives (see also eq. |B9| ) 



Ik 



I fec| 3 

4 |fcc| 3 



2D 
3D 



(13) 



We now compare the expectations from linear response 
theory to the numerical simulations. We begin by tak- 
ing the Fourier transform of SB from our 2D numerical 
simulations from a downstream region where the shock 




kc /" P 

Fig. 3. — Spectral evolution of magnetic field from the slice at 
840c/o;p from Fig. [2] Initial field spectrum (black solid line) is 
plotted after 450 uip 1 (red), 900 uip 1 (green), and 1350 ujp 1 (blue) 
based on simulation data. Dashed curves represent analytic evolu- 
tion of the initial field and overpredict decay of short-wavelength 
structures. 



is fully developed, i.e., behind the shock front at x = xo, 
1 



SB(x ,y) = 



/2w 



dk y SB ky (x ) exp (ik y y) (14) 



We evolve SBk y in accordance to equation (fT2")) using the 
asymptotic forms in equation (|13[) and compare our ana- 
lytically evolved spectra with the numerical simulation 
at later times. In Figure [3l we take the initial data 
from a region at xo = SAQc/uo p , which we marked with 
a line in Figure [5] To accumulate sufficient statistics, 
we average the fields over 14 c/u> p in the flow direction. 
We evolve these spectra for 450 (red), 900 (green), and 
1350 LJp 1 (blue) using equation (fT3|) and compare this 
to snapshots taken from our numerical simulations at 
these times. Theory and simulation agree at very low 
wavenumber (k y c/u) p < 0.2). However, theory overpre- 
dicts the cutoff in power at larger k. The discrepancy 
suggests that linear theory is insufficient to describe the 
nature of downstream magnetic turbulence and that ad- 
ditional physics is needed (see 0. 

The lack of power in short wavelength modes suggests 
that the total B-field is determined by long wavelength 
modes. We now use equation (fl"2"|) to find a simple decay 
law for the total B-field: 



5B(x ,y,t) 



1 



dk y SB ky (x ) exp (ik y y - j k t) , 



Q 5 ) 

where SBk (x) is defined in equation (|14[) . Inserting 
equation (fT3"|) into (fl5)) and integrating, with SBk y = ak'P, 
where p is the long wavelength spectral index and a nor- 
malizes the amplitude, we find 



SB(x ,y,t)-- 



1 



a I dkyk^exp 

27T Jk 



-aui„t 



kyC 




aujpt / 



j(16) 

(p+1)/3 

I (17) 
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Fig. 4. — Magnetic energy density (in units of upstream kinetic 
energy) as a function of position downstream of the shock. A 
broken power law proportional (x — aipeak) -2 ^ 3 tits well at early 
times, but a (x — x pe ak) — 1 power law fits better at later times. 



where we have taken y = without loss of general- 
ity, a — 4/-7T in 3D and a — 1 (2D), fco is small such 
that oj p t(koc/u!p) 3 ~ 0, and T is the gamma function. 
Note that lu p is for the downstream plasma frequency. 
However, this detail is irrelevant in terms of determin- 
ing the power law as a function of t. Setting ko to be 
small is a safe approximation in the simulations, and 
of course is excellent in much larger astrophysical sys- 
tems. Thus if the initial spatial spectrum is a power law 
in wave number |5-Bfc| 2 ex k 2p , then the theory predicts 
SB 2 cx t- 2 (P+!)/ 3 . 

We study the region immediately following the peak 
of magnetic energy in the shock front (where it reaches 
SB 2 = i? 2 lax ) at £poak and plot its value as a function 
of position in the postshock region. For a shock mov- 
ing at constant velocity, we have a; pe ak — x oc t. Hence 
SB 2 oc (zpoak — x)~ 2 ( p+1 ^ 3 . Our numerical simulations 
are extremely suggestive that the magnetic energy den- 
sity follows the p = 0, e B = SB 2 /Sir cx t~ 2 / 3 decay 
expected for an initially flat magnetic spectrum at early 
times, then steepens to a t~ x decay at later times as 
shown in Figure [U until the noise finally swamps the sig- 
nal. 8 This t^ 1 decay rate would imply an initial spatial 
index of p — 1/2 at low wavenumber, though our analysis 
of additional simulations with a large transverse spatial 
scale suggest p = 0. The difference in the index of the 
decay law expected from the theory and measured from 
the simulations may be due to magnetic trapping (see 
Sj4j. Gruzinov (2001b) offered an alternative explanation 
of t^ 1 decay of the magnetic energy density observed in 
2D simulations of Weibel instability in counterpropagat- 
ing pair plasmas. In that simulation the beams were 
moving perpendicular to the simulation plane 9 . In or- 
der to explain the field decay, Gruzinov (2001b) had to 

8 Gruzinov (2001a) performed shock simulations (initiated by a 
collision between two e plasmas), which show decay of the mag- 
netic energy averaged over the dimension across the flow similar to 
the early phases of the decay we report here. 

9 Being orthogonal to the direction of motion these simulations 
did not form a shock and retained some counterstreaming at late 
times. In contrast, our simulations lose the counterstreaming once 
the shock forms. 



assume that the magnetic structures increase their size 
at the Alfven velocity. In our shock simulations, the 
downstream magnetic structures do not significantly ex- 
pand (§2), hence the similarity of decay law between the 
two simulations is likely a coincidence. The picture of 
filament expansion and merging of Gruzinov (2001b) is 
likely valid very close to the shock, but is not observed 
in the longer evolution of the downstream plasma. 

4. MAGNETIC TRAPPING 

Our simulations and theory suggest power law t~~ q with 
q ~ 1 temporal decay of the total magnetic energy den- 
sity in rough agreement with each other. However, linear 
theory and numerical simulations disagree on the decay 
rate for short wavelengths, which hinders a determina- 
tion of the ultimate fate of the fields on times longer 
than several thousand plasma periods. This discrep- 
ancy may arise from the nonlinear effects of magnetic 
trapping. Particles in the magnetic fields do not follow 
straight line trajectories that are weakly perturbed, but 
are partially trapped and strongly deflected. We illus- 
trate this point from following test particles' orbits in 
numerical simulations in Figure [5l The test particles 
suffer large deflections from straight-line orbits as they 
encounter magnetic clumps. Therefore, the Larmor radii 
of many of the test particles are of the same order of the 
sizes of these clumps or smaller. Indeed closer inspec- 
tion of some of the test particle orbits that are embed- 
ded inside the clumps suggests that they are completely 
trapped. 

These strong departures from weakly perturbed parti- 
cle dynamics may be the cause of the decreased damping 
at large wavenumber found in the simulations, compared 
to the predictions of unmagnetized plasma theory. Mag- 
netic trapping may also modify the decay of magnetic 
fields at small wavenumber, leading perhaps to a different 
decay law instead of expected t 2 / 3 decay law that we find 
from our simple linear theory (eq. |17| with p = 0). The 
temporary (and permanent, for some particles) binding 
of particles to the spatially intermittent magnetic fields 
reduces the effect of rapid phase mixing central to the 
unmagnetized plasma damping theory. Magnetic trap- 
ping already plays a critical role in the saturation of the 
initial Weibel instability at the shock transition region 
(Kato 2005; Davidson et al. 1972). Our analysis suggests 
that trapping also plays an important role downstream 
even though the average magnetic amplitudes are greatly 
reduced from the shock transition region - the essential 
point is that within the isolated filaments, the magnetic 
pressure is not small. The problem of the damping of iso- 
lated magnetic structures in a plasma with partial mag- 
netization illustrated in Figure [5] will be the subject of a 
separate investigation. 

5. DISCUSSION 

We have studied the downstream evolution of mag- 
netic turbulence in the context of a collisionless e + e~ 
shock both analytically and numerically. Our large scale 
2D simulations show the formation of filaments in the 
foreshock region which merge and grow until they reach 
the shock transition region. Past the shock transition 
region, these filaments break up into magnetic clumps 
in a quasi-homogenous medium where the background 
particle distribution function is an isotropic Maxwellian. 
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Fig. 5. — Segments of particle orbits from the high resolution 2D simulation overplotted on a snapshot of magnetic energy downstream 
of the shock. Typical sizes of magnetic clumps are 10-20 c/u> p . The shock is near the right boundary of this slice, moving to the right. The 
large angle deflections in particle orbits suggest the Larmor radii of many particles in this slice are of the same order as the sizes of the 
magnetic clumps. 



In such a background, we showed that magnetic energy 
will decay like t~ q with q ~ 1 due to untrapped parti- 
cle phase mixing, which is broadly consistent with our 
numerical simulations. In detail, the theoretical decay 
rates depend more strongly on wavelength than is seen 
in the numerical experiments. We suggested that mag- 
netic trapping may play an important role in resolving 
this discrepancy. Trapping can lead to MHD-like be- 
havior and possible long time persistence of some of the 
magnetic energy in spatially intermittent, more strongly 
magnetized subregions. This effect is a subject of further 
study. 

Magnetic held decay may be alleviated via an inverse 
cascade from small scales to large scales. An inverse cas- 
cade via current filament merging operates strongly in 
the foreshock region (Silva et al. 2003; Spitkovsky 2005). 
However, it is unclear if this picture can be applied to the 
downstream region. The magnetic filaments so prevalent 
in the upstream region are gone and are replaced by iso- 
lated magnetic clumps (loops), a result common to both 
2D and 3D simulations. Filament merging does occur, 
but is confined to the foreshock, where the filaments ex- 
ist. We see little evidence of merging magnetic clumps. 
This is reinforced by Katz et al. (2007), who suggest via 
self-similar arguments that the inverse cascade does not 
operate downstream of the shock. Hence, the magnetic 
energy should damp away in the downstream region. 

We have not studied in detail how these magnetic 
clumps are confined in the downstream plasma. How- 
ever, Fujita et al. (2006) have found these same mag- 
netic clumps in simulation of the non-relativistic Weibel 
instability that is appropriate for galaxy clusters. In this 
case, they argue that pressure of the external medium 
confines these magnetic clumps. The strong perturba- 
tion from straight line motion, which we studied in Sj2]for 
test particles, suggest that the magnetic clumps in our 
simulation are also sensitive to momentum transfer be- 
tween particles and itself. That is, the magnetic clumps 
in our simulations are also confined by external pressure. 

Whatever the final fate of the spatially intermittent 
fields, our study of the nature of e + e _ shocks suggests 
that the belief in the persistence of Weibel generated 
magnetic fields with strengths comparable to those ap- 
pearing within the shock transition is overly optimistic. 
Magnetic field energies tend to decline rapidly after 



about a few hundred plasma skin depths. Only the very 
long term evolution is still open to some question, i.e., 
does (es) settle at some value smaller than 10~ 3 , or does 
the spatially intermittent magnetic field decline to zero? 

The rapid decay suggested by our analysis puts severe 
constraints on the synchrotron emission mechanisms for 
GRBs. The width of the emitting region is ~ 10 9 c/w p 
(Piran 2005b), which is much larger than the region over 
which we expect magnetic fields to persist. However, 
decaying magnetic fields may not be inconsistent with 
GRB observations. Pe'er & Zhang (2006) suggest that 
the prompt emission from the internal shock may be 
more consistent with a decaying magnetic field compo- 
nent rather than a persistent field component. A field 
that persists over a scale of 10 4 — 10 5 plasma skin depths 
fits the spectra better at low-energies than a field that 
persists over the entire thickness of the shell (10 9 plasma 
skin depths). Small magnetized regions may also be im- 
portant in the context of the afterglow (Rossi & Rees 
2003). 

In young pulsar wind nebulae, post shock magnetic 
fields averaged over the whole latitudinal extent of the 
observed emission tori have energy densities of a few per- 
cent of the post shock plasma energy density. It is possi- 
ble that weaker magnetic fields exist near the midplanc 
of the equatorial flow, a region of particular interest to 
the conversion of flow energy into the observed nonthcr- 
mally emitting spectra of e . If so, the Weibel mediated 
shock dynamics studied here may be of relevance to these 
systems' behavior. 

It is also possible that weak systematic upstream mag- 
netic fields are of essential importance, and that shocks in 
completely unmagnetized plasmas are an oversimplifica- 
tion. Suprathermal particle generated at the relativistic 
shock front may alter the basic physics of the collision- 
less shock, if a mean field is present. Milosavljevic and 
Nakar (2006b) argue accelerated particles streaming into 
the upstream magnetized medium can drive long wave- 
length, magnetized turbulence with SB/B 3> 1, which 
might persist into the downstream and provide the mag- 
netization required in phcnomcnological models of GRB 
and PWN emission. Then the shock mediated by Weibel 
turbulence becomes a subshock within a much larger ex- 
tended structure, responsible only for thermalizing the 
bulk of the flow and injecting the particles that are Fermi 
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accelerated in the turbulence generated by high energy 
particle streaming. This is a relativistic version of Bell's 
(1978) (also see Bell 2004, 2005) picture of particle ac- 
celeration in non-relativistic shocks. 

Finally, we briefly mention the extension of this theory 
to electron-ion plasmas. Let us presume initially that 
the ions and electrons are isotropic but remain decou- 
pled, i.e., a two-temperature relativistic plasma. Equa- 
tion (fTTjl becomes 



X 



■ 7T u p,i + ^p,e 

'4 \k\cLU ' 



(18) 



where = Amnie 2 /jiirii is the plasma frequency of 
the ions, 7$ is the Lorentz factor associated with the 
ion temperature, to 2 e = Aitn e e 2 /j e m e is the plasma fre- 
quency of the electrons, and 7 e is the Lorentz factor as- 
sociated with the electron temperature. Depending on 
the relative values of the electron temperature and ion 
temperature, one term may dominate. However, initial 
large-scale simulations of ion-electron collisionless shocks 
suggest that both reach roughly equipartition with each 
other (Spitkovsky 2007), thereby reproducing the physics 
of the shock. In this case, the relativistic electrons 



and ions contribute equally to the decay rate because 
thermal equipartition prevails, i.e., m e ^ e — rriiji. Thus, 
the electron-ion plasma has the same dynamics as the 
electron-positron plasma. 



We thank S. Cowley, D. Kocelski, M. Milosalavjic, 
A. Pe'er and E. Quataert for useful discussions. P.C. 
and J. A. thank the Institute for Advanced Study for 
its hospitality; P.C. also thanks the Canadian Institute 
for Theoretical Astrophysics for similar hospitality. P.C. 
is supported by the Miller Institute for Basic Research. 
J. A. has benefited from the support of NSF grant AST- 
0507813, NASA grant NNG06GI08G, and DOE grant 
DE-FC02-06ER41453, all at UC Berkeley; by the De- 
partment of Energy contract to the Stanford Linear Ac- 
celerator Center no. DE-AC3-76SF00515; and by the 
taxpayers of California. A.S. is pleased to acknowledge 
that the simulations reported on in this paper were sub- 
stantially performed at the TIGRESS high performance 
computer center at Princeton University which is jointly 
supported by the Princeton Institute for Computational 
Science and Engineering and the Princeton University 
Office of Information Technology. 



APPENDIX 

SUSCEPTIBILITY IN TWO AND THREE DIMENSIONS 



In this section, we solve the susceptibility (eq.[9]) for two and three dimensional plasma. The 3D case ( §A.1|) has 
been previously solved by Mikhailovskii (1979) and we reproduce his result here for completeness. In addition, we 
present the 2D case ( ^A.2[l . which allows for a comparison to the two-dimensional simulations reported in this paper. 



Three Dimensions 
After summing over the electrons and positrons, equation §§§ is 

d kv x 
,dp x 



4?rx 



"p,NR 



fa 



d 3 P , 



(Al) 



kVy dp y t 

where n is the number density of electrons and positions. Since /o is isotropic and therefore independent of angle, we 
orient the spherical integral in a non-standard manner so that the pole points along the k- vector, i.e., the y-axis. We 
find p ■ y = cos 9 and p ■ x = sin 9 sin cf>. So equation (|A1|) becomes 

,2 



4ttx 



v sm sm' 



1 + 



$0 2 1 jo 

— — p dpdil, 
dp 



where dfl = sin 9d9d<f>- Performing the integral over 



Lu/kv — cos 9 
and making the substitution ( 



cos 9, we find 



4ttx 



7T0J, 



p,NR 



v(l-( 2 ) 1 + 



c 



uj/kv — £ 



Integrating over £ from -1 to 1, we find 

4ttx 



2jrW p,NR 




p dp 



1 / l+uj/kv 

2 ° S \l-u/kv 



U3 

kv 



dfo 2 7 ■> 
— p dpdQ. 
dp 



(A2) 



(A3) 



(A4) 



Note that the logarithmic function will give an imaginary part when cu r /kcf3 < 1. This implies that waves whose phase 
velocity, oj r /k, is small compared to the thermal speed of the background particles, cf3, will be damped (also see - 
We pull this imaginary component out of the equation, which makes this damping more explicit, to find: 



4ttx = 



2 ™p,NR 



dfo 



p dp—v{ — - — 



dp 



kv I 



kv 



1 



f-V 

\kv) 



1, fl+iv/kv 
- log — 

4 & \l-u)/kv 



i-e (k 2 v 2 
2 v 



(A5) 
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where 8 is the unit step function (i.e., Q(z) = 1 for 5i(z) > and Q(z) — for 3t(z) < 0). Equation (| A5|) is precisely 
Mikhailovskii (1979)'s result. We now apply a three dimens iona l relativistic Maxwellian fo oc exp (—E/kT), which is 
appropriately normalized, J d 3 pfo — n, and solve equation (|A5|) numerically. For the ultrarelativistic case v sa c, we 
find a simple form in the limit u r <C kc: 

^"'iw- (A6) 

4 \k\cuj 
Two Dimensions 

Starting from equation (jAlj) . we assume fo = f(p2d)g(Pz), where p2d = \Jpx + Py an d perform the integral over p z . 
The resulting two dimensional analogue of equation (lAljl is 

. w p,NR f ( d kv x d \ / 2 /a-n 

4ttx = ^- / TT + - — r— 3— -dp, (A7) 



UJ — kVy dp y J Tl 

where we have dropped the subscript "2c?" from p. Defining p x — pcosO, p y — psm9, and similarly for v x and v y , we 
find: 

4 7TX = ^m [paipf* devco S >e(i+ , , (As) 



(uj/kv) — sin 9 J dp 

We may transform the 9 integral from to 2n to a contour integral over the unit circle by making the appropriate 
substitutions (Carrier, Krook, & Pearson 1983) 

cos6>^ i (z + z" 1 ) ,sin0-> (z - z" 1 ) (A9) 

dJd^,r^L (A10) 
^ Jo Jr 

where T is the unit circle. After a bit of algebra, we find 



The singular points in this equation are z± = (iuj/kv) =p ^/l — (uj/kv) 2 . We perform the contour integral by noting 
that only the Z- root contributes for (cu/kv) 2 < 1: 



: 2m 4NR /"./, ("\ 2 df 



4nx = _;_£^ / V 1 " ^ ^PdP- (A12) 
wfcn J V \ kv/ dp 

We apply a two dimensional relativistic Maxwellian / oc exp (— E/kT), where / is appropriately normalized, i.e. 
J fd 2 p = n. Expanding to lowest order in to/kv, we find: 

uj 2 

A ^~ l TU^- ( Al3 ) 
\k\cuj 

DECAY OF AN INITIAL FIELD OF FLUCTUATIONS 

The plasma susceptibilities (cq. 1 1 1 j ) define the linear response of the plasma. These susceptibilities are complex and 
hence the electric permittivity e = 1 +47r% is also complex. The implications of Poynting's theorem for the propagation 
of small amplitude waves and fluctuations in such a medium have been described in many texts and monographs (e.g., 
Bekefi 1966, Melrose and McPhedran 1991, Stix 1992). For completeness, we give a brief discussion of the theory 
behind expression (|12p . 

We begin with Poynting's theorem: 



dt V 8tt 

where S = (c/Att)E x B is the Poynting vector. We assume field energy is distributed uniformly in our uniform 
plasma, so V • S = 0, i.e., there is no transport of energy spatially. In addition, for nonpropagating Weibel modes, 
6E <C SB. Thus, we find a simplified expression: 

Writing the fields with truncated amplitudes 

dV T v{r, t)-< Qt g ( _ T/2j r/2)j | r j ^ ^ (Ui) 
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with SBtv — when the coordinates are outside of the volume, V, and time is outside of the interval (— T/2, T/2). 
Thus, we can define space-time averages of the fields while still expressing them in terms of convergent Fourier 
transforms. We write the Fourier transforms as 

fijku — d 3 k dudjTV exp (ik ■ r — iujt) . (B4) 

J J — oo 

We express the field energy density and the Sj ■ SE work in terms of the Fourier amplitudes and then average over 
time and space. Taking T and V to oo, we integrate over the resulting (5-functions to find 

§1 (^) = / d 3k J d»\{&3to6EL + c -c->, (B5) 

where () represents the space time average over the fluctuation wavelengths and variability times and c.c. is the complex 
conjugate. We apply the same Fourier transform to 8B and apply the same average over time and space. We use the 
linear response (eq.(S]) 8E ku — (i/x^)Sjku to find 



1 d(\5B ku) \ 2 ) 
8tt dt 

Now using 5jkui = ikcSBkuj/^n, we find 

d{\5B kuJ \ 2 ) 



dt 

where jk w is the damping rate or 



-S^x)" 1 ^! 2 ). (B6) 



-2lU\SB ku) \ 2 ), (B7) 



(/,c)a -9f(47r X )- 1 . (B8) 



Equation (|B7|) is a specific form of the fluctuation-dissipation theorem (Thompson and Hubbard 1960) from which the 
same result can be derived. 

Equations (|B7p and ()B8|) define the evolution of magnetic energy in terms of the linear response. For notational 
simplicity we reduce ku> — > in the rest of the text. We obtain simple forms for 7/., using the asymptotic forms of Anx 
from equation (fTTj) for 3D and 2D. We find for k<w p : 

^=<3*3D- (B9) 



The cubic dependence on wavenumber in equation (|B9[) suggests short wavelength modes are very strongly damped, 
while longer wavelength modes can survive much longer. 

Taking the general form of iirx for 2D and 3D from equation (|A5|) and (|Al2j) , we numerically compute the damping 
rate from equation (IB8j) for Weibel modes where cu r = 0, because of the non- propagating nature of the downstream 
magnetic chimps. We show the results in Figure [B6l Also plotted are the simple forms for 7^ from equation (|B9[) . 
Though, the asymptotic forms are only valid for kc/uj p -C 1, the numerical and asymptotic results are in excellent 
agreement throughout. Thus for simplicity, we use equation (|B9[) (also eq.[T3] in the main body) for the damping 
rates. 
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l)A5|l and l)A12Jl_ for 3D (thick solid line) and 2D (thin solid line) respectively. Also overplotted are the analytic forms of these expressions 
from equation H13H for 2D (thin dashed line) and 3D (thick dashed line). 
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